Pathways of bond topology transitions at the interface of silicon nanocrystals and 

amorphous silica matrix 
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The interface chemistry of silicon nanocrystals (NCs) embedded in amorphous oxide matrix is 
studied through molecular dynamics simulations with the chemical environment described by the 
reactive force field model. Our results indicate that the Si NC-oxide interface is more involved than 
the previously proposed schemes which were based on solely simple bridge or double bonds. We 
identify different types of three-coordinated oxygen complexes, previously not noted. The abundance 
and the charge distribution of each oxygen complex is determined as a function of the NC size as 
well as the transitions among them. The oxidation at the surface of NC induces tensile strain 
to Si-Si bonds which become significant only around the interface, while the inner core remains 
unstrained. Unlike many earlier reports on the interface structure, we do not observe any double 
bonds. Furthermore, our simulations and analysis reveal that the interface bond topology evolves 
among different oxygen bridges through these three-coordinated oxygen complexes. 

PACS numbers: 61.46.Hk, 68.35.Ct 
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I. INTRODUCTION 

After a long arduous effort, photolumincsccncc from 
silicon has been achieved from its nanocrystalline 
form A critical debate, however, continues over the 
nature of the interface chemistry of silicon nanocrys- 
tals (Si-NCs) embedded in amorphous silica which has 
direct implications on the optical activity of the inter- 
face @, M, EL 0, H- Wolkin et al. reported that the ox- 
idation of porous silicon quantum dots results in a red 
shift in the photolumincsccncc (PL) spectra which in- 
dicates the importance of oxygen-related interface bond 
toplogy [2]. Along this line, Puzder and co-workers com- 
pared PL calculations of nanoclusters with different pas- 
sivants and surface configurations and proposed the main 
reason for the red shift to be double Si=0 bonds 
Countering this, Luppi et al. reported excitonic lumi- 
nescence features caused by Si-O-Si bridge bonds at the 
surface of silicon nanoclusters 1- As a supporting evi- 
dence for the latter, Gatti et al. recently demonstrated 
that Si-O-Si is the most stable isomer configuration ||. 
To reconcile, Vasiliev et al. claimed that bridge bonds 
and double bonds have similar effect on PL Q. 

All of the work cited above represent density functional 
theory (DFT)-based calculations with small Si clusters of 
less than 100 atoms surrounded by either passivants like 
hydrogen or oxygen [2, 3] . But actual samples are 

profoundly different: the fabricated systems consist of Si- 
NC with a diameter larger than 1 nm, embedded in amor- 
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phous silica (a-SiO^) matrix. Identifying this fact, Tu 
and Tersoff studied Si/a-Si02 interface using a bond or- 
der dependent empirical potential and proposed that Si- 
O-Si bridges are the main blocks at the interface lowering 
the surface energy [7j. Using the same model potential 
Hadjisavvas et al. studied Si-NCs embedded in a-Si0 2 
and they also reported bridge bonds as the mechanism 
for lowering surface strain energy Q. In our previous 
work on the formation and structure of Si NCs, we also 
observed that this model potential falls short to charac- 
terize the system accurately j_especially the structure and 
chemistry at the interface [9(. As a rigorous approach 
to interface structure and dynamics, Pasquarello et al. 
used first-principles molecular dynamics (MP) technique 
to investigate planar Si/Si02 interface [1(|. They ob- 
served that the oxygen atoms momentarily are bonded 
to three silicon atoms during oxidation process. The net 
effect of these threefold coordinated oxygen atoms dur- 
ing oxidation process was to expel silicon atoms out of 
interface. This is interpreted as a balancing process to 
decrease the increased density of Si/Si02 interface due 
to oxidation However, no substantial breakthrough 
was made over the last decade on this issue. 



II. METHOD 

In this work, we employ the reactive force field 
(ReaxFF) model developed by van Duin et al. which 
improves Brenner's reactive bond order model [Tl| to a 
level of accuracy and validity allowing molecular dynam- 
ics simulations of the full reaction pathways in bulk [l2| . 
The parameters for this force field were obtained from 
fitting to the results of ab initio calculations on rele- 
vant species as well as periodic boundary condition DFT- 
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based calculations of various crystalline polymorphs of 
relevant materials. The ReaxFF calculates bond orders 
which is the measure of bond strength from local geom- 
etry. This allows realistic chemical environment such as 
over/under coordination and bond breaking/formation 
for large-scale (about 5000 atoms) MD simulations. 

To facilitate our discussion regarding the surface- 
bonded oxygen complexes, we distinguish among three 
different types of silicon atoms. We label those silicon 
atoms with all silicon neighbors each with zero oxidation 
state as c, to denote core silicon atom. Among the re- 
maining (non-c) silicon atoms, those with at least one 
bond to c are labeled as s, denoting as a surface sili- 
con atom. For further investigation of NC, we seperate 
core Si atoms into two sub-categories as inner-core and 
outer-core atoms: Among core Si atoms which have at 
least one surface Si neighbour categorized as outer-core 
Si atoms and rest of core Si atoms categorized as inner- 
core Si atoms. Finally, any other silicon atom is labeled 
as m, denoting matrix silicon atom. Hence, a complex 
consisting of an oxygen atom bonded with two surface sil- 
icon atoms is labeled as ss. The other oxygen complexes 
are sm, ssm, sss, mms as sketched in Fig. [Tj where the 
last three correspond to three-coordinated oxygen (3cO) 
atoms. 

We use ReaxFF to represent the interactions in the 
model system. We start with a large simulation cell (box 
length 43 A) of silica glass formed through a melting and 
quenching process used by one of us [H, [3| earlier to 
study silica glasses earlier. Next, similar to Hadjisavvas 
et al. Q, we delete all atoms within a predetermined 
radius to insert crystalline silicon to form NC. In this way, 
we create NCs with radii ranging from 5.5 A to 16.7 A . 
For the largest NC we insert 967 Si atoms into a spherical 
hole with radius 16.7 A created in amorphous matrix. 
Even for this case, the minimum distance between NC 
surface to simulation box face is about 5 A which can 
still accommodate the interface layer. We also pay special 
attention in the removal of spherical region so that the 
correct stoichiometry for the amorphous matrix is met. 
Thus, our amorphous matrix has two O atoms for every 
Si atom with a density of 2.17 g/cm 3 which is the density 
of glass at room temperature and atmospheric pressure. 

We set periodic boundary conditions in all directions 
and while keeping NC at 100 K, we employ simulated an- 
nealing process to Si02 region to end up with an amor- 
phous matrix free of artificial strain around the NC. Then 
we set whole systems' temperature to room temperature 
(300 K) and continue performing MD simulation for 75 ps 
to have thermal equilibrium between the two regions. We 
set MD simulation time step to 0.25 fs for all simulations. 
For every 62.5 fs time interval, we have record the config- 
urations to analyze the transitions taking place between 
different bond topologies. The numbers attached to each 
arrow in Fig. [T] indicate the total number of registered 
transitions during the simulation in that direction be- 
tween the complexes for a representative NC of radius 
13.4 A. Almost balanced rates in opposite directions is 
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FIG. 1: (Color online) The transitions between different oxy- 
gen complexes bonded to the interface. Dark green (dark 
gray) large spheres represent matrix silicon atoms, and the 
light green (light gray) large spheres represents surface sili- 
con atoms of the NC, small red (dark gray) spheres represent 
oxygen atoms. The numbers indicate the number of transi- 
tions recorded in the simulation in each direction among the 
complexes for the NC of radius 13.4 A. 

an assurance of the attainment of the steady state in 
our simulation. Note that we do not observe any direct 
transition other than the paths indicated in Fig. [TJ For 
instance, a direct transition of the complex ss to sm does 
not take place, but it is possible through an intermediate 
transition over the ssm which is a 3cO. We should also re- 
mark that the balanced transitions continue to take place 
after the steady state is attained which indicates that the 
interface bond topology is dynamics, i.e. not frozen. 

To analyze Si-NC/a-Si02 interface, we construct NC 
surface using Delaunay triangulation scheme (Fig. [2] in- 
set) [15| . In two dimensions a Delaunay triangulation of 
a set of points corresponds to defining triangles such that 
no point in the set is inside the circumcircle of any trian- 
gle. In three dimensions triangles extend to tetrahedra 
and circumcircles become circumspheres. Since our NCs 
are nearly spherical in shape, we triangulate projection 
points of surface Si atoms onto the unit sphere. Hence, we 
apply Delaunay triangulation over the two dimensional 
6-<f) plane. We can then create surfaces using this triangu- 
lation. This surface enables us to calculate every atoms' 
distance to NC surface. By this means we plot various 
data such as charge, bond order etc. with respect to dis- 
tance to the surface, to extract information about surface 
chemistry of Si-NC embedded in amorphous matrix. 

III. RESULTS AND DISCUSSIONS 

For different NC radii we observe similar trends in 
bond order distribution, average charges etc., therefore, 
we present only the figures of the system for a typical 
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FIG. 2: (Color online) Radial distribution function of Si 
atoms in NC. Inset: First peak resembles Si-Si bond length 
distribution centered around 2.34 A with 0.3 A FWHM value. 
Second peak's width resembles bond angle deviation. Inset: 
Representation of NC surface created with Delaunay Trian- 
gulation. Blue dots represents surface silicon atoms. 



NC of radius 13.4 A. A useful data to elucidate the struc- 
ture of these systems is the radial distribution function 
(RDF). In Fig. rj we present RDF of NC atoms only, 
where the first broad peak centered around 2.34 A with 
a 0.3 A full width at half maximum value (FWHM), rep- 
resents Si-Si bond length distribution in NC (Fig. [2] in- 
set). The maximum extent of the NC can also be read 
from same plot at about 27 A, where the RDF goes to 
zero. Observation of a broad first peak at Fig [2] demands 
further investigation of Si-Si RDF of NC atoms. For 
this purpose, in Fig. [3] we present bond length proba- 
bility distributions (akin to RDF) for three categories of 
Si atoms: inner core (with no bonds to surface atoms), 
outer core (bonded to surface) and surface NC atoms. 
We observe from Fig. [3] that Si-Si bond lengths in the 
inner core are centered around the equilibrium value and 
have a narrow width due mainly to thermal vibrations, 
whereas the bond length distributions of outer core and 
surface atoms have increasing shift for the most proba- 
ble bond length values and broader widths. These shifts 
and particularly the increase in distribution widths can- 
not be attributed to thermal broadening. Taken together 
these two observations is a clear indication of increasing 
strain as a function of distance from the center of the NC. 
To further investigate this deviation of Si-Si bonds from 
crystalline Si, in Fig. [5] we present bond length distribu- 
tion with respect to distance to NC surface averaged over 
2 ps of simulation time after the steady state is reached. 
This figure illustrates the gradual development of radial 
strain from the center to NC surface. These observations 
show clearly that oxidation at the surface of NC results in 
a tensile strain at Si-Si bonds which becomes significant 
only around the interface, while keeping the inner core 
almost unstrained. This tensile strain in the NC agrees 
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FIG. 3: (Color online) Bond distance probability distribution 
of Si atoms in NC. Solid line represents inner core Si-Si bonds, 
dashed line outer core Si-Si bonds, and dotted line represents 
the surface Si-Si bonds. 

with previous measurement of Hofmeister et al. 16]. 

Another consequence of this tensile strain is that the 
total bond orders of core-NC atoms are somewhat smaller 
than those of oxide-Si's as seen in Fig. In the same fig- 
ure we also show the calculated net charges using electron 
equilibration method [13] • Nearly zero net charges of the 
core-Si atoms reflects the covalent type of bonding well 
within the NC. The bonding becomes increasingly ionic 
away from the NC core as observed by the charges of Si 
atoms which reach the value of 1.3e at the oxide region 
(Fig. [5]). As a result, the positive charges of surface-Si 
atoms form a shell at the surface of NC. This observa- 
tion is similar with those obtained with DFT calculations 
[l8l [T^]. On the other hand, negative charges of oxygen 
atoms bonded to surface form another shell that enclose 
NC and finally total average charges approach to zero 
within the oxide region. In Fig. [5] we also observe that 
the magnitude of average charges of oxygen atoms which 
are bonded to surface are greater than those in the ma- 
trix. But, the bond orders are nearly the same. This is 
due to existence of 3cO atoms bonded to surface. Note 
that the average bond order of oxygen atoms which are 
bonded to surface is about two (cf., Fig. [5]). Thus, those 
oxygen atoms form three partial bonds, two strong and 
one weak bond. Finally, we would like to note that un- 
like many others 0, [H, IB] , we do not observe any double 
bonds. 

The occurrence of 3cO has been noted by a number of 
groups. Pasquarello showed that the bistable E[ defect of 
ev-quartz structure may lead to 3cO as a metastable state 
as well as Si-Si dimer bond and calculated the energy of 
the former to be higher than the latter. Pasquarello pro- 
posed that 3cO acts as an intermediate metastable state 
during structural relaxations at the interface [2fj| . Simi- 
larly Boero et al. observed 3cO atoms in their ab initio 
calculations [2l[ and reported this feature as a metastable 
state. In Table U we present the collected statistical data 
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TABLE I: Statistical results of atom charges and numbers, N, for all NC diameters, -Dj^q considered. Abbreviations for atom 
types are explained in Fig. [T] Charges are in the units of electronic charge and the angles 8, are in degrees. 
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FIG. 4: (Color online) Variation of Si-Si bond length averages 
(calculated over 1 A wide bins) as a function of distance from 
the NC surface -which is defined by Delaunay tesselation. The 
solid line is a fit to the data to guide the eye. 
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FIG. 5: (Color online) Top: Silicon bond orders (triangles) 
and charges (circles) as a function of distance from surface of 
NC; Middle: Oxygen bond orders and charges as a function 
of distance from the surface of NC; Bottom: Total average 
charge as a function of distance from the surface of NC. The 
averaging bin width is 1 A. 



at the end of the simulation of 75 ps. For all oxygen com- 
plexes, the number of bridges, average charges of bridge 
oxygens, and the average bridge angles for s-O-s are tab- 
ulated. We observe in Table U that the number of sss 
complexes is very small due to narrow bond angle re- 
quirement of this configuration. For the remaining 3cO 
complexes, ssm and mms, their percentages are seen to 
increase with curvature. This can explain the fact that 
other studies [20L l2lT ] which have identified the 3cO com- 
plexes as metastable were based on the planar Si/SiC>2 
interfaces. So, this is an indication of the importance of 
curvature in the stability of 3cO complexes. Hence, as 
one would expect, there is a linear relation between the 
total number of bridges with surface area as indicated in 
Fig. [SJ This finding is supported by Kroll et al. who re- 
ported 3 and 33 such bridges for Si-NC with radii 4.0 A 
and 7.0 A, respectively [l8j ] . 



IV. CONCLUSIONS 

In conclusion, the realistic chemical environment pro- 
vided by reactive force field model enables us to under- 
stand the bond topology of Si-NC/a-Si02 interface and 
its internal dynamics. Particularly, it reveals that there 
are different types of oxygen complexes at the Si-NC sur- 
face some of which contain 3cO complexes whereas there 
are no double bonds. The curvature has a positive effect 
on the occurence of 3cO. The relative abundance of dif- 
ferent complexes and their charge and geometrical char- 
acteristics are extracted. The inner core is observed to be 
almost unstrained while the outer core and the interface 
region of the NC are under increasing strain up to about 
a few per cents. In general, our work clearly shows that 
the Si NC-oxide interface is more complicated than the 
previously proposed schemes which were based on solely 
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FIG. 6: (Color online) The number of bridges at the Si-NC 
surface vs radius squared. The line is a linear fit to data. 



simple bridge and double bonds. The provided informa- 
tion here paves the way to construct realistic Monte Carlo 
moves for the simulation of large-scale silicon nanostruc- 
tures embedded in oxide matrix. 
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